Optimized Effective Potential Method in Current-Spin Density Functional Theory 
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Current-spin density functional theory (CSDFT) provides a framework to describe interacting 
many-electron systems in a magnetic field which couples to both spin- and orbital-degrees of free- 
dom. Unlike in usual (spin-) density functional theory, approximations to the exchange-correlation 
energy based on the model of the uniform electron gas face problems in practical applications. In this 
work, explicitly orbital-dependent functionals are used and a generalization of the Optimized Effec- 
tive Potential (OEP) method to the CSDFT framework is presented. A simplifying approximation to 
the resulting integral equations for the exchange-correlation potentials is suggested. A detailed anal- 
ysis of these equations is carried out for the case of open-shell atoms and numerical results are given 
using the exact-exchange energy functional. For zero external magnetic field, a small systematic low- 
ering of the total energy for current-carrying states is observed due to the inclusion of the current in 
the Kohn-Sham scheme. For states without current, CSDFT results coincide with those of spin density 
functional theory. 
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I. INTRODUCTION 

Density functional theory (DFT) [1, 2] is the standard 
tool to calculate the electronic structure of interacting 
many-electron systems. The original theorems of DFT 
have been successively extended to cover a wide range 
of different physical situations. One of the most im- 
portant extensions is the spin-DFT (SDFT) formalism [3] 
which allows to describe systems with magnetic ground 
states in arbitrary external magnetic fields. In SDFT, 
however, the magnetic field orJy couples to the elec- 
tronic spin while the coupling to the orbital degrees of 
freedom is not taken into accoimt. To include also the 
coupling to the orbital currents, one has to resort to the 
current-spin-density functional theory (CSDFT) of Vig- 
nale and Rasolt [4, 5]. 

Conceptually, DFT, SDFT, and CSDFT are very simi- 
lar: they all map the system of interacting electrons onto 
a system of non-interacting particles in some effective 
fields. In the case of DFT, this auxiliary system yields 
the same electron density as the interacting one while in 
SDFT the magnetization densities of the two systems co- 
incide as well. In CSDFT, also the paramagnetic current 
density of the auxiliary system is equal to the paramag- 
netic current density of the true, interacting system. In 
all three formalisms the energy of the interacting system 
is written as a functional of the corresponding densities 
and the value for the ground state energy is obtained by 
minimizing this functional with respect to the densities. 

While all three flavors of DFT are exact in principle, 
in practice they all require an approximation for the 
exchange-correlation (xc) energy (which is a piece of the 
total energy) as a functional of the respective densities. 
Both in DFT and in SDFT, approximations based on the 
uniform electron gas such as the local (spin-) density 
aprroximation (L(S)DA) are surprisingly successful and 
they also form a good starting point for the construction 
of more sophisticated fimctionals such as generalized- 
gradient approximations (GGA's). 



In CSDFT, however, the situation is different. While 
it is possible to construct LDA-type approximations 
along similar ideas as in (S)DFT, these functionals are 
awkward to use in practical calculations [6] for a clear 
physical reason: when a uniform electron gas is ex- 
posed to an external magnetic field. Landau levels form 
and, for given magnetic field, the xc energy density ex- 
hibits derivative discontinuities as function of the den- 
sity whenever a new Landau level starts to be filled. If 
this xc energy density is then used as the main ingre- 
dient in the construction of an LDA, these discontio- 
nuities show up in the corresponding xc potentials at 
those points in space where the local densities coincide 
with the densities of the uniform gas for which these 
discontinuities occur. This makes practical caluclations 
extremely difficult. One solution to this problem is to 
smoothly interpolate the xc energy density between the 
limits of weak and strong magnetic fields [7, 8] but this 
interpolation then misses the physics in the exchange- 
correlation energy arising from the Landau levels. 

The problem described above is entirely due to the 
construction of the LDA in CSDFT. Since the appearance 
of Landau levels is intrinsically an orbital effect, the use 
of explicitly orbital-dependent approximations to the xc 
energy functional offers a promising alternative which 
we explore in the present work. In (S)DFT, orbital func- 
tionals have attracted increasing interest in recent years 
[9] since these approximations offer a cure for notori- 
ous problems like, e.g., the self-interaction error of lo- 
cal and semilocal functionals. The calculation of the xc 
potential corresponding to an orbital-dependent xc en- 
ergy functional is technically somewhat more involved 
than for explicitly density-dependent approximations 
and can be achieved with the so-called Optimized Ef- 
fective Potential (OEP) method [10]. Here we present 
the extension of the OEP method to the CSDFT formal- 
ism [11] and derive a set of coupled OEP integral equa- 
tions for the corresponding exchange-correlation poten- 
tials. We then introduce a simplifying approximation in 
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the spirit of Krieger, Li, and lafrate (KLI) [12, 13] which 
transforms the integral equations into a set of algebraic 
equations which can be solved more easily in practical 
calculations. We present the resulting equations for fully 
non-coUinear effective magnetic fields. This generalizes 
earlier work on a non-collinear implementation of or- 
bital functionals in SDFT [14]. It is also similar in spirit 
to a recent work [15] which, however, uses a much larger 
set of densities. 

In a previous paper [11] we have solved the OEP 
equations of CSDFT for two-dimensional quantum dots 
in external magnetic fields. In the present work we 
study the case of atoms at zero external magnetic field 
using the exact-exchange energy functional. In partic- 
ular, we are interested in studying open-shell atoms 
which generally have degenerate ground states. It is 
well-known [16, 17] that in (S)DFT common approxi- 
mations for the exchange-correlation energy do not lead 
to the same total energy for the different ground-state 
configurations. Here we investigate this problem in the 
framework of CSDFT because some of these ground 
states actually carry a non-zero current and one might 
hope that CSDFT is better suited to describe the de- 
generacies. The motivation to study simple atomic sys- 
tems also aims at a better understanding of the differ- 
ences between SDFT and CSDFT. The same orbital func- 
tional may perform quite differently in the three differ- 
ent schemes. 

The present paper is organized as follows: in Section 
II we give the fundamental equations of CSDFT before 
we present the derivation of the non-collinear OEP- and 
KLI-equations in Section III. This is followed in Section 
IV by an analysis of the properties of the KLI-potentials 
for open-shell atoms at zero external magnetic field. Nu- 
merical results are presented in Section V before we 
draw oui conclusions. 



II. CURRENT-SPIN DENSITY FUNCTIONAL THEORY 

In this Section we briefly describe the formalism of 
current-spin density fimctional theory as originally sug- 
gested by Vignale and Rasolt [4, 5]. With this extension 
of the original DFT [1, 2] formulation it becomes pos- 
sible to study interacting many-electron systems in ex- 
ternal magnetic fields. The CSDFT approach also goes 
beyond the widely used spin-DFT formalism [3] in the 
sense that while in SDFT the magnetic field only couples 
to the spin degrees of freedom, in CSDFT it also couples 
to the orbital degrees of freedom through the vector po- 
tential. 

The Hamiltonian of interacting electrons in an exter- 
nal electrostatic potential vo{r) and an external magnetic 
field Bo(r) = V x Ao(r) is given by (atomic units are 



used throughout) 

H = f + W + Jd^r h{r)vo{r) - Jd^r ih{r)Bo{r) 

-t-i ydVjVWAoW-h^^ ydVn(r)Ag(r),(l) 

where T and W are the operators of the kinetic energy 
and the electron-electron interaction, respectively. The 
operators for the density, the magnetization density, and 
the paramagnetic current density, are given by 



n(r) = ^'1'(r)*(r) , 
m(r) = -/XB^^(r)<7*(r) 



(2) 
(3) 



and 



jp(r) = ^ (i'\r)V^{r) - (V^'t(r))«'(r)) , (4) 

respectively. Here we have defined field operators 
^^(r) = (i/)|(r), V'|(r)) for two-component spinors, i.e., 
the formulation is not restricted to coUinear magnetism 
with all the spins aligned in a single direction, a is the 
vector of Pauli matrices and /is = l/(2c) is the Bohr 
magneton. 

Following Vignale and Rasolt [4, 5], the ground state 
energy can be written as a functional of the three densi- 
ties as 

E[n, m, jp] = F[n, m, jp] -|- j d^r vo{r)n{r) 
- Jd\ m(r)Bo(r) + i Jd^r jp{r)Ao{r) 
^ ld'rn{r)Al{r), (5) 



+ 



where F[n, m, jp] is a imiversal fimctional of the densi- 
ties n, m, and jp in the sense that it is independent of the 
external fields vq, Bq, and Aq. It may be decomposed in 
the usual way as 

i^[n, mjp] = Ts[n,in,jp] + U[n] + E^c[n,m,jp] , (6) 

where Tg[n,in,jp] is the kinetic energy functional for 
non-interacting electrons. 



U[n] = i jd'r jd' 



, n(r)n(r') 



(7) 



is the classical electrostatic energy and Exc[n,m,^j^ is 
the exchange-correlation energy. 

Using Eq. (6) and minimizing the total energy (5) 
leads to the effective single-particle Kohn-Sham (KS) 
equations of CSDFT which read 



ei$.(r), (8) 
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where $i(r) are two-component, single-particle Pauli 
spinors. The effective potentials are given by 

Vs{r) = voir) + VH{r) + v,c{r) + ^ [aI{v) - Al{r)] , 

(9) 



and the paramagnetic current density by 



B«(r)=Bo(r)+B^e(r) 



and 



As(r) = Ao(r) + Aj;c(r) , 
where the Hartree potential is given by 

n(r') 



VH{r) = j^^r' 



r — r' 



(10) 



(11) 



(12) 



The exchange-correlation potentials are functional deri- 
vatives of the exchange-correlation energy E^c with re- 
spect to the corresponding conjugate densities 



Bxc(r) 



5Exc[n,ni,jp] 
5n{r) 



6E, 



[n,m,jp\ 



Aa;c(r) = c 



(5m(r) 



SExc[n,m,jp 



(13) 



(14) 



(15) 



The effective potentials are such that the ground-state 
densities of the Kohn-Sham system reproduce those of 
the interacting system. The particle density can then be 
computed by 



N 



(16) 



i=l 



the magnetization density by 



N 



m(r) = -/XB ^i>|(r)cr$j(r) 



(17) 



i=l 



N 



W = ^ E H WV$i(r) - (V$l(r)) $,(r)] , (18) 



where the sums in Eqs. (16) - (18) run over the occupied 
Kohn-Sham spinor orbitals. 

As usual in DFT, the exact form of the exchange- 
correlation energy functional Exc[n,m.,jp] is unknown 
and needs to be approximated in practice. In the present 
work, we focus on a class of approximate function- 
al which also has attracted increasing interest in SDFT 
in recent years. This is the class of functionals which 
explicitly depend on the Kohn-Sham orbitals and are 
therefore only implicit fimctionals of the densities. In 
our context these functionals are appealing for two rea- 
sons: first, they are constructed without any reference 
to the uniform electron gas and second, they are ide- 
ally suited to describe the appearance of Landau levels 
which in itself may be viewed as an orbital effect. 
III. OPTIMIZED EFFECTIVE POTENTIAL METHOD IN 
CSDFT 

The calculation of the exchange-correlation potential 
for orbital-dependent functionals in ordinary (S)DFT is 
done in the framework of the Optimized Effective Po- 
tential (OEP) method [9, 10, 18]. This method derives 
its name from the fact that it yields that local potential 
whose orbitals minimize the total energy fimctional of 
the interacting system. This optimized potential is ob- 
tained as a solution of the so-called OEP integral equa- 
tion. 

Recently, the OEP equations have been derived for the 
non-collinear formulation of SDFT [14]. Here we derive 
the OEP integral equations for the exchange-correlation 
potentials in CSDFT [11] by calculating the functional 
derivatives of E^c with respect to the three effective po- 
tentials Vs, Bs, and A.,. Making use of the correspon- 
dence both between Kohn-Sham spinors and ground- 
state densities, as well as between Kohn-Sham spinors 
and effective potentials, these functional derivatives can 
be computed in two different ways by using the chain 
rule, i.e.. 



5Exc 



6Vs{t) C OVs{v) OVs{v) 



N 



5Exc S^ijr') 



+ h.c. 



(19) 



5Exc 



^B^r) + c^"^^'^ ^5B,(r) ^^^^"^ ^ 5B«(r) 



N 
i=l 



6Exc S<P,{r') 
5M^') ^B.(r) 



h.c 



(20) 



SEx 



SAs{r) c SAs{r) SAs{r) 



N 



6$i{r') SA,{r) 



h.c 



(21) 
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Eqs. (19) to (21) constitute a system of coupled integral 
equations for the unknown exchange-correlation poten- 
tials. For simplicity, we have assumed that the approxi- 
mation for to be used depends only on the occupied 
spinor orbitals such as, e.g., the exact-exchange func- 
tional 



E: 



EXX 



^ occ 



1,3 



(22) 



which we use in our numerical implementation. 

For any approximation of E^c given explicitly in 
terms of the spinor orbitals, the functional derivatives of 
Exc with respect to these spinors can be evaluated easily. 
The other functional derivatives in Eqs. (19) - (21) may 
be computed exactly from first-order perturbation the- 
ory by considering variations of the Kohn-Sham spinors 
due to a perturbation SHs{r') of the Kohn-Sham Hamil- 
tonian. To first order in the perturbation these variations 
are 

(5$i(r) = V / dV'$](r')5i?«(r')$i(r') , (23) 

M» I 



where for simplicity we have assumed a non-degenerate 
spectrum. For arbitrary variations 5vs{v), 5^s{^), and 
5 As (r) in the three effective potentials, the perturbation 
5Hs (r) is given by 



5Hs(v) = 5vs{y) + —W6As{y) 
Zic 

+ ^5A,(r)V + ^As{v)5As{v) + iiBcr5^s{v) .(24) 

Insertion into Eq. (23) allows us to identify the func- 
tional derivatives of the spinors with respect to the ef- 
fective potentials as 



E 

3 = 1 



$](r')$i(r') 



(25) 



5$i(r) 



JB,(r') 



(26) 



3 = 1 



and 



6As{v') 



E [$](r')V'$,(r') - (V'<l>](r')) $,(r')] + ^A,(r)$t(r')<l>.(r')| • (27) 



From Eqs. (25) - (27) one can compute the static response 
functions, i.e., the functional derivatives of the densities 
with respect to the effective potentials. Inserting every- 
thing into Eqs. (19) - (21) one can then write the OEP 
integral equations in a very compact form as 



N 



^$t(r)*,(r) + /i.c. = 0, 



N 



-liB $l(r)o-*i(r) + h.c. = 



(28) 



(29) 



and 



1 ^ 

- {$l(r)V*i(r) - [V$l(r)] *,(r)} + h.c. = , 

i=l 

(30) 

where we have defined the so-called orbital shifts [9, 19] 



3 = 1 



(31) 



with 



+ i-A,e(r) [$](r')V'$i(r') - (v'<i>](r')) $.(r') 



2ic 

t) E 

+MsB,e(r')$kr')o-$i(r') - ^^r'^ — 



<5$l(r') 



(32) 



The name "orbital shifts" (31) derives from their struc- 
ture as a first-order shift from the unperturbed Kohn- 
Sham orbital <l>i under a perturbation whose matrix el- 
ements are given by Dij . The OEP equations (28)-(30) 
have a very simple interpretation: they merely say that 
the densities do not change xmder this perturbation. 
Keeping in mind that the Kohn-Sham system already 
yields the exact densities, this statement is actually quite 
obvious. 

As already mentioned, the OEP equations are a set of 
coupled integral equations for the exchange-correlation 
potentials. In this work we do not attempt a full solution 
of these equations but rather suggest a simplifying ap- 
proximation [18] in the spirit of Krieger, Li, and lafrate 
(KLI) [12, 13] who introduced the same approximation 
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in the usual OEP method of SDFT. In this KLI approxi- 
mation the orbital shifts are approximated by replacing 
the energy denominator by some constant, i.e., 

*'(^) i - $i(r)Ai j . (33) 

Inserting this approximation into the OEP equations 
and applying the completeness relation for the Kohn- 
Sham spinors one obtains after some algebra a set of 
algebraic equations for the exchange-correlation poten- 
tials which can conveniently be written as 

P(r)Vxc(r) = n{v) . (34) 

Here, we have defined the 7-component vector Vxc(r) as 

^--^ (35) 



V4(r)= (^^.c(r),B^,(r),-A^,(r) 

and the 7x7 matrix V{r) has the structure 

/ n(r) -m'^(r) jj(r) \ 
V{v) = -m(r) /x|n(r)l J{r) , 
\ j,(r) J^{r) M{r) J 

where 1 is the 3x3 unit matrix. The matrix elements of 
the 3x3 matrices J and TV are defined by 



Ja/3(r) 

and 



1=1 \ 



(36) 



1 ( 



dvcc drp drp dra 
1 dn{r) dn{r) 



4n(r) dra dr/s 



(37) 



where a = 1, 2. 3 corresponds to the cartesian coordi- 
nates X, y, and z, respectively. The seven components 
of the vector Tl{r) on the right-hand side of Eq. (34) are 
given by 

T ^ f SE \ 

1 ^ 

^i+a(r) = ^E 

^ i=l 

( - fiB'^l{r)aa4^7 + mi,a{r)Dii + h.c] ,(39) 
1 d<^l{r) SE,. 



with the density rii (r), magnetization density (r), and 
paramagnetic current density jp ^(r) of the single orbital 
$i(r). It is worth mentioning that in order to arrive at 
this result we used the identity [5] 



(41) 



V n(r)A,e(r) =0 



which follows directly from gauge invariance of the 
exchange-correlation energy. 

The KLI equations (34) can be solved by iteration: 
start with an intial guess for the potentials to compute 
the orbitals and the right-hand side of Eq. (34), then 
solve this equation for the new potentials and iterate im- 
til self-consistency is achieved. 

In a different work [11] we have solved the KLI equa- 
tions for a two-dimensional quantum dot in an external 
magnetic field. In the present work we use our CSDFT- 
OEP formalism to study open-shell atoms in zero exter- 
nal magnetic field. In the next Section we discuss some 
further assumptions we employed in our implementa- 
tion and deduce some analytic results for the KLI poten- 
tials. 



IV. OPEN-SHELL ATOMS AT ZERO MAGNETIC FIELD 

We want to employ our CSDFT-OEP formalism to 

study open-shell atoms. From the point of view of CS- 
DFT this is interesting since some states out of the mul- 
tiplet of degenerate ground states have a non-vanishing 
current density while others do not carry a current. 

In the limit of zero external magnetic field, the Kohn- 
Sham equation (8) takes the form 



— + voir) + VH{r) + v,c{r) + ^ (VA,e(r)) 



+ -A^eV + HB<TB^c{r) *i(r) = ei$i(r). (42) 
ic I 



We further employ the coUinear approximation assum- 
ing that the Kohn-Sham spinors decompose into spin- 
up (cr = +1) and spin-down (a — —1) orbitals, i.e., 
$,(r) = (^,,,=+i(r),0) or $(r) = (0, (^(r),,,=_i). As 
a result the magnetization density is parallel to the z- 
direction, m(r) = (0, 0,m(r)). In addition, we assume 
cylindrical symmetry for both the densities and the cor- 
responding conjugate potentials (that is they do not de- 
pend on the azimuthal angle (/)). As a consequence the 
magnetic quantum number is a good quantum number 
for the single-particle orbitals which then take the form 



(r) = Qiair, 6) exp{im(j))x{(j) 



(43) 



2i dra S<^l{r) 



+ 3p,kA^)Dkk + h.c. \ (40) 



where we used radial coordinates and m is the mag- 
netic quantum number (not to be confused with m(r), 
the z-component of the magnetization density), a is 
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the spin index and x(cr) is the eigenfunction of the z- 
component of the spin operator. In the coUinear approx- 
imation, Bj;c(r) = (0, 0, Bxc(jy) is parallel to the z-axis 
while Aa:c(r) = ^xciy)^<^ where is the unit vector in 
(^-direction. As an additional consequence of the cylin- 
drical symmetry of our problem we have VA ,;c(r) = 0. 

We restrict ourselves to ground states whose densities 
can be reproduced by a single Slater determinant. For 
example, for the boron atom one configuration has all 
three up-electrons and the two down-electrons in states 
with magnetic quantum number to = while in another 
configuration one of the up-electrons occupies an m = 1 
state with the other occupations imchanged. In this way 
current-carrying and zero-current states can be consid- 
ered. The resulting current only has a component in 
the (/)-direction, jp(r) = jp{r)e^. We may then rewrite 
Eq. (42) as 



1 TO 

— -I- Vo{r) + VH{r) + Vxcir) H — A^c(r) 

2 cr sm t) 



(r) = (r).(44) 



In the following we discuss a number of typical cases: 
For atomic closed-shell configurations, where the den- 
sity is spherical and both the magnetization density 
and the paramagnetic current density vanish identically, 
both Axe and B^c vanish identically. Obviously, in this 
situation CSDFT reduces to the original density-only 
DFT. 

For groimd-state configurations where only orbitals 
with 171 = are occupied, the correct value for jp{r) 
- which is zero at any point in space - is trivially ob- 
tained already within the SDFT scheme. Therefore we 
expect that ^^.^(r) = ^^(r), B^dr) = B^^^^{r), and 
^xc(r) = 0. Actually, any other choice of Axc{i^) would 
not make any difference for the groxmd state densities. 
In a way this may be regarded as a simple manifestation 
of the non-uniqueness of the CSDFT potentials pointed 
outinRef. [20]. 

As a third case we consider ground-state configura- 
tions with a half -filled shell as in, e.g., the nitrogen atom. 
Again, SDFT already gives the correct values of the to- 
tal densities. Therefore, we again expect that Vxd'r) = 
v'^r^^ir), B^^v) = B^J^^ir) and AM = 0. 

Ground-state configurations carrying a non-vanish- 
ing paramagnetic current are the most interesting ones 
in our context. At zero external magnetic field, this 
situation only arises for open-shell atoms away from 
half-filling. Indeed, it is for these states that we expect 
Axcir) ^ as well as v^dr) ^ t^ff ^^(0, ^.^(r) ^ 

In the following we analyze the KLI equations for 
the above cases in order to confirm these expectations. 

We remind the reader that in our derivation of the OFF 
equations we assumed that E^c depends only on the oc- 



cupied orbitals. Moreover, we also assume that 



' exp{—imq 



(45) 



holds. Both of these assumptions are correct for the 
exact-exchange functional which we consider in our nu- 
merical implementation. 

Under these assumptions, the first two KLI-equations 
(for a = ±1) are 

VxcA^) + -^^^.e(r) = + 4?,.(r) , (46) 

c ncr(r) 

where we have defined the spin-dependent scalar po- 
tential 

Vxc,a{r) = Vxc{r) + i^B(rBxc{r) ■ (47) 
The terms on the right-hand side of Eq. (46) are given by 



^^laW = T^'^nima{r)Uxc,imair) (48) 



and 



w. 



(2) 



with 



and 



1 SEx 



(49) 



(50) 



(r) (Vxc,a{r) - Uxczmcrir)) 



+ ^ jp,ima{r)Axc{r) . (51) 

Here nimcr(r) and jpAma{Y) are the density and the para- 
magnetic current density of the single orbital ipima{^) 
which, for our symmetry, are related by 



r smf 



(52) 



The third KLI equation reads 

{^jpA^>xcA^)^ + \N{r)Axc{v) = 

Y.{w%{v)+w'il{v)) (53) 



with 



occ -2 ( \ 
. 'Hma\*- ) 



(54) 
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and 



(55) 



(56) 



It is interesting to note that Vxc,air) and A^dr) couple 
to each other orvly if at least one of the jp,cr(r) is non- 
vanishing. 

At this point, we again consider open-shell config- 
urations for which all occupied orbitals have m = 0. 
Then iV(r) vanishes and the KLI equation (46) reduces 
to the KLI equation of SDFT while Eq. (53) becomes a 
trivial identity. As a consequence, ^-^^..^(r) = ?;,f,^/'-^(r) 
and j4j;c(r) is undetermined. As discussed above, AxdT^) 
does not affect any of the ground state densities and we 
fbc it as Axc{r) = 0. 

Next we consider configurations with a half-filled 
shell. We assume that we have already solved the SDFT 
KLI equations and use the resulting orbitals and poten- 
tials plus the initial guess Axdr) = as a start for the 
iterative solution of the CSDFT KLI equations. We sub- 
stitute this initial guess into Eqs. (46) and (53) to com- 
pute the new potentials. The occupied orbitals of SDFT 
either have m = or they come in pairs with m and — m. 
This leads to the same contributions to Uxc.a,i{^) for or- 
bitals in the same shell while for the paramagnetic cur- 
rent they contribute with equal magnitude but opposite 
sign. Hence, the KLI equations become 



w 



(1) 



' (r) 



W 



(2) 



Ur)=v!,^rir) {57) 



and 



:7V(r)^-"'(r) = 0^A:r(r)=0 



(58) 



This shows that the SDFT solution along with Axc{y) = 
is also a CSDFT solution. We also tested numerically 
that the solution A^d^) = is stable against (not neces- 
sarily small) perturbations of the initial guess. 

Finally, we consider the most interesting case of 
groimd-state configurations with a paramagnetic cur- 
rent different from zero. For these configurations we 
expect Ar^d^) 7^ 0. Solution of the third KLI-equation 
(53) with respect to A^d^) yields 



Axd^) = 

E(7 (w'xc,a(r) + W)icla(r) - ip,a(r)Wxc,a(r) 

' W) 



.(59) 



In this equation, the denominator increases for increas- 
ing number of electrons. The numerator also typically 
increases when more orbitals are occupied but, due to 
large cancellations for contributions arising from or- 
bitals with opposite values of m, it increases with a 



slower rate than the denominator. As a consequence, 
we expect larger exchange-correlation vector potentials 
^xc(r) for atoms in the first row than for atoms in the 
second row (but the same column) of the periodic table. 

In the remainder of this Section we discuss the asymp- 
totic behavior of the exchange-correlation potentials and 
the vector potential. 

We start by assuming that, for finite systems, the 
exchange-correlation potentials in the asymptotic region 
far away from the system behave as 



and 



lim ^xc(r) = 



(60) 



(61) 



Eq. (60) certainly is a reasonable assumption: it is 
the weU-known asymptotic behavior for v^ca of SDFT 
which we assume to be xmchanged in CSDFT. Eq. (61) 
then ensures that the term proportional to A^d^)/''' if* 
the Kohn-Sham equation (44) decays faster than W:rc,cr(r) 
asymptotically. At this stage, Eq. (61) may be viewed 
as a working assumption in order to be able to proceed 
further. Below we show that it is consistent with the so- 
lution of the KLI equation. 

Under this assumption we can deduce [9, 21] the 
asymptotic behavior of the atomic orbitals from the 
Kohn-Sham equation (44) as 



lim ¥'im(7(r) 



-i+i/ft, 



(62) 



where f3im(7 = \J—'^Uma- This implies that the spin den- 
sity is dominated asymptotically by the highest occu- 
pied orbital of that spin. The same is true for the current 
density if the magnetic quantum number of this orbital 
is different from zero (as typically is the case for current- 
carrying states of open-shell atoms). 

In order to proceed further with our analysis we 
restrict ourselves to the exact-exchange functional of 
Eq. (22). Then the KLI equation (46) allows us to estab- 
lish the following relation between v^ca and A^^ in the 
asymptotic region 

lim ?^xc,c7(r) + --^^r^^xc(r) ---l-rfxc,Ar„M„CT , (63) 
r^oo c r sm u r 

where we tacitly assumed that we are taking the limit 
away from a nodal plane of the highest occupied or- 
bital of spin (7 [9, 22]. Here No- is the orbital index 
of that orbital and is its magnetic quantum num- 
ber. Since we are working in the coUinear approxima- 
tion, the Kohn-Sham equations for the two spin chan- 
nels are completely decoupled and we can choose a con- 
stant shift in Vxc,a such that 







(64) 



Eq. (63) together with Eq. (60) then implies 



which is consistent with the assumption of Eq. (61). 

However, a closer inspection of the KLl equations 
(46) and (53) shows that they become linearly depen- 
dent in the asymptotic region and therefore do not have 
a unique solution. This again may be viewed as a 
manifestation of the non-imiqueness problem in CSDFT 
[20, 23]. In our numerical procedure to be described 
in the next Section we take a pragmatic approach to 
the problem of linearly dependent KLl equations and 
choose a solution with ^xc(r) and a Vxcai'"^) satisfy- 
ing Eq. (60). 

Before concluding this Section, we discuss some sym- 
metry properties of the exchange-correlation vector po- 
tential and exchange-correlation magnetic field. By in- 
spection of the two KLl-equations (46) and (53) it is clear 
that imder the exchange of spin-up and spin-down elec- 
trons, Bxc(r) changes sign. Similarly, exchanging an 
electron from an orbital with magnetic quantum num- 
ber TO to a previously unoccupied one with —to leads to 
AxciTc) changing sign. These transformations can be per- 
formed independently leading to the same total energy. 



V. NUMERICAL RESULTS FOR OPEN-SHELL ATOMS 
AT ZERO EXTERNAL MAGNETIC FIELD 

In this Section we describe the numerical results for 
open-shell atoms obtained within the KLl approxima- 
tion of CSDFT using the exact-exchange functional of 
Eq. (22). In particular, we are interested in calculating 
total energies of current-carrying and zero-current states 
in various ground-state configurations. In principle, the 
states of the ground-state multiplet should be degener- 
ate but in SDFT zero-current states are always lower in 
energy than current-carrying states. Since the current 
appears to be the quantity leading to these spurious en- 
ergy splittings, it is interesting to see if CSDFT (where 
the current is one of the fundamental variables) can al- 
leviate the problem. Since the main difference between 
SDFT and CSDFT is the appearance of an exchange- 
correlation vector potential in the Kohn-Sham equation, 
we also present some results for the xc potentials in the 
different approaches. 

We have developed an atomic code for CSDFT and 
SDFT calculations in a basis set representation, assum- 
ing cylindrical symmetry of the Kohn-Sham potentials 
and densities. As basis functions we use Slater-type 
functions for the radial part multiplied with spherical 
harmonics for the angular part. For the Slater func- 
tions we employ the quadruple zeta basis sets (QZ4F) 
of Ref. [24]. We have tested our code by computing the 
total energies of spherically symmetric atoms of the first 
and second row of the periodic table in exchange-only 
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FIG. 1: Spherical component of the exchange vector potentials 
for current-carrying states of the oxygen and sulfur atoms. 



KLl approximation and compared with results from ac- 
curate, fully numerical codes available in the literature 
[9, 12, 25]. Our code reproduces these total energies 
to within an average deviation of 0.1 kcal/mol for the 
atoms in the first row and to within an average devia- 
tion of 0.5 kcal/ mol for atoms in the second row of the 
periodic table. 

We performed self-consistent exchange-only calcula- 
tions in the KLl approximation of CSDFT for open- 
shell atoms in current-carr5dng and zero-current con- 
figurations. The configurations are selected by speci- 
fying the number of occupied states for each value of 
the magnetic quantum number to. Once a choice has 
been made, the occupations remain unchanged during 
the self-consistency cycle. In all the cases we studied 
we were able to obtain self-consistent solutions for both 
zero-current and current-carrying states. 

As expected, for zero-current states we always ob- 
tain a self-consistent CSDFT solution with vanishing ex- 
change vector potential, (r) = 0. In fact, this solution, 
which is equivalent to the corresponding SDFT solution, 
always gives the lowest total energy. 

For current-carrying states we always find a non- 
vanishing Ax- However, as a consequence of the linear 
dependence of the KLl equations (46) and (53) discussed 
in the previous Section, the exchange potential and vec- 
tor potential are not uniquely determined. In fact, with- 
out additional numerical measures we obtain rmphysi- 
cal potentials which diverge asymptotically. This may 
lead to a wrong ordering of the occupied and unoc- 
cupied orbitals or even to convergence problems. A 
nirmerically convenient scheme to use the KLl equa- 
tions (46) and (53) during the self-consistency cycle is 
to first calculate t'a:c,<7(r) from Eq. (46) with the A^dr) 
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FIG. 2: Spherical component of the exchange scalar potentials 
for the current-carrying state of the oxygen atom computed in 
SOFT and CSDFT. 



obtained in the previous iteration. In this step we im- 
pose the asymptotic behavior of i^rccrlr) ^^^^ —1. Then, 
with this updated Vxcai^) we use Eq. (53) to obtain a 
new AxcX^^). In order to enforce the asymptotic hmit 
^a;c(r) '^^f we add a small quantity 5 to iV(r) in Eq. 
(59). Total energies and current densities are very insen- 
sitive to the choice of 6: for the fixed value oi S = 10""' 
a.u. which we use for all our calculations, total energies 
vary by an order of 10^^ kcal/mol or less if S is varied 
by an order of magnitude aroxmd its chosen value. 

In Fig. 1 we show the L = (i.e., spherical) com- 
ponent of an expansion of the exchange vector poten- 
tials in terms of Legendre polynomials, i.e., Ax{r) = 
X^^^Q^^(r)Pi:,(cos^), for the oxygen and sulfur atoms 
in the current-carrying state. As we argued in the pre- 
vious Section, the exchange vector potential of sulfur is 
smaller than the one for oxygen which also implies that 
the difference between SDFT and CSDFT total energies 
is smaller for the heavier atom. 

In previous work [26-28] it has been assumed that 
CSDFT and SDFT calculations lead to very similar re- 
sults because the term associated with the vector po- 
tential is expected to be small. In order to verify this 
assumption, we compare the self-consistent exchange 
potentials and magnetic fields for the current-carrying 
state of the oxygen atom for a SDFT and a CSDFT cal- 
culation. The spherical components of the exchange po- 
tentials in the two approaches are shown in Fig. 2. The 
potentials are hardly distinguishable on the scale of the 
plot which confirms the initial assumption. If one looks 
at the corresponding exchange magnetic fields (Fig. 3) 
one sees that the overall structure of the SDFT and the 
CSDFT results are quite similar. However, there are sig- 
nificant differences in magnitude close to the nucleus 
which can be expected to have a visible effect on the re- 
sulting chemical shifts [29]. This is also reflected by a 
substantial difference in the relative magnetization den- 
sity C(0) = (nT(0) - ni{0))/{n^{0) + 71^(0)) at the nu- 



clear position. For the current-carrying state of oxygen 
we obtain the value ((0) = -1.03 x IQ-^ in SDFT and 
C(0) = -1.16 X 10-3 in CSDFT which amounts to a dif- 
ference of approximately 13%. While we are confident 
that the difference of SDFT and CSDFT values for ^(0) 
is not a numerical artefact, the absolute numbers should 
be read with some caution. In order to estimate the accu- 
racy of these numerically sensitive results we have also 
calculated the same quantity for the nitrogen atom and 
obtain ^(0) = —1.77x10-^. This value differs by approx- 
imately 9% from the value C(0) = -1.62 x 10"^ given 
in Table 10 of Ref. [9] which was obtained with a fully 
numerical code for spherically symmetric effective po- 
tentials. 

Finally, we have calculated the spurious energy split- 
tings between different configurations of open-shell 
atoms of the first two rows of the periodic table in SDFT 
and CSDFT. The results are collected in Table I. As 
mentioned before, in all cases the splittings are posi- 
tive, i.e., the zero-current states are always lowest in en- 
ergy. The splittings themselves are systematically lower 
in CSDFT than in SDFT due to the additional varia- 
tional degree of freedom in the former approach. The 
effect of including the vector potential is more signifi- 
cant for the lighter atoms while for the second row the 
CSDFT splittings are only less than 0.1 kcal/ mol smaller 
than the ones obtained from SDFT. Although a CSDFT 
approach to the degeneracy problem appeared promis- 
ing our results show only a small and insufficient im- 
provement. This finding is somewhat at odds with re- 
cent suggestions to reduce the splittings by inclusion of 
the current density as a variable in the construction of 
approximate exchange-correlation functionals [30-32]. 
These works, however, suggested orbital functionals in 
the framework of SDFT where the orbital dependence 
entered through the current density. No attempt was 
made to implement these functionals in a CSDFT frame- 
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FIG. 3: Spherical component of the exchange magnetic field 
for the current-carrying state of the oxygen atom computed in 
SDFT and CSDFT. 
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Atom 


^SDFT 


^CSDFT 


Atom 


^SDFT 


y\CSDFT 


B 


1.66 


1.38 


Al 


1.68 


1.58 


C 


1.58 


1.34 


Si 


1.76 


1.63 


O 


2.36 


2.29 


S 


3.04 


3.01 


F 


2.32 


2.27 


CI 


3.15 


3.10 



TABLE I: Spurious energy splittings A = E{M = ±1) — 
E{M = 0) (in kcal/mol) between current-carrying and zero- 
current states computed in SDFT and CSDFT. 

work. Finally, we point out that in a recent work [33] 
we showed, that a pure DFT approach using the exact- 
exchange functional surprisingly leads to energy split- 
tings which are almost an order of magnitude smaller 
than the ones obtained in SDFT. 



VI. CONCLUSIONS 

In this work, we have shown how one can use orbital- 
dependent functionals in the framework of current-spin 

density functional theory. We have derived the OFF 
integral equations which have to be solved to obtain 
the corresponding exchange-correlation potentials. We 
have simplified these integral equations in the spirit of 
the well-known KLl approximation. 



We have analyzed the KLl equations and the result- 
ing potentials for open-shell atoms at zero external mag- 
netic field and have also presented numerical results for 
these systems using the exact-exchange functional. We 
have shown that CSDFT and SDFT are equivalent for the 
states with zero paramagnetic current. This equivalence 
breaks down for current-carrying states where total en- 
ergies in CSDFT are lower than those of SDFT. 

We also verified that the CSDFT Kohn-Sham scheme 
leads to a reduction (compared to SDFT) of the spurious 
splittings between current-carrying and zero-current 
states although it is too small to recover the degeneracy 
between these states. 

The most important result of our study, however, 
is the fact that the problems of LDA-type current- 
density functionals derived from the uniform electron 
gas (such as tinphysical discontinuities of the corre- 
sponding exchange-correlation potentials) never appear 
when using orbital dependent functionals. 
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